package scu.maqiang.possion;
import java.util.Arrays;

import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T61;
import scu.maqiang.mesh.Mesh2T6;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;
public class Possion2DT61 {

		public static void main(String[] args) {
	        int N = 8, label = 10;
	        double ep = 1.0 / N;
	        
	        //ʹ��lambda���ʽ�����ί��
	        ScalarFunc regionFunc = (xy, llabel, param) -> {
	            double xCell = xy[0] * N - (int)(xy[0] * N);
	            double yCell = xy[1] * N - (int)(xy[1] * N);
	            return (Math.abs(xCell - 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25)? 1.0:-1.0;
	        };

			Mesh2T6 mesh2D = new Mesh2T6().square(32, 32, T3Type.Left);
			//mesh2D.displayNodes();
			//mesh2D.displayElements();
			//mesh2D.displayBoundarys();
	        mesh2D.setDomainLabel(regionFunc, null, label);
	        
	        FES2T61 fs = new FES2T61(mesh2D);
	        SRMatrix A = new SRMatrix(fs.GetNdof());

	        ScalarFunc regionCoef = (xyz, relabel, param) -> relabel == label ? 0.01 : 1.0;
	        ScalarFunc ConstCoef = (xyz, relabel, param) -> 1.0e0;
	        ScalarFunc[] gcf = {regionCoef};
	        gcf[0] = ConstCoef;
	        
	        ScalarFunc coefFunc = (xyz, relabel, param) -> 2 * xyz[0] * (1 - xyz[0]) + 2 * xyz[1] * (1 - xyz[1]);
	        fs.assembleStiff(gcf, null, BVPType.COMMON, A);
	        //fs.assembleHeatStiff(new double[] {1.0}, BVPType.COMMON, A);
//	        System.out.println(A);
//	        System.out.println(A.makeEntrySymmetric());
	        double[] RHS = new double[fs.GetNdof()];
	        //GeneralCoefFunc ConstCoef = (xyz, relabel, param) -> 1.0;
	        gcf[0] = coefFunc;
	        fs.assembleSource(gcf, null, BVPType.COMMON, RHS);
	        A.makeEntrySymmetric();
	        fs.applyBCZero(A, RHS, Direct.X, BCType.MBN, 1, 2, 3, 4);
	        
//	        System.out.println(A);

//	        IterSSolver solver = new IterSSolver(1.0e-15, 1000000);
//	        
	        double[] x = new double[fs.GetNdof()];
//	        CSRMatrix mat = new CSRMatrix(A);
//	        solver.PCGSSOR(A, RHS, x, 1.5e0, 1);
//	        Arrays.fill(x, 0.0);
//	        solver.PCGSSOR_2(A, RHS, x, 1.5e0, 1);
	//
//	        
//	        //System.out.println(A);
//	        
//	        //System.out.println(mat);
//	        
//	        Arrays.fill(x, 0.0);
//	        solver.CG(mat, RHS, x, 1);
//	        
//	        Arrays.fill(x, 0.0);
//	        solver.CG(A, RHS, x, 1);
	        
	        Arrays.fill(x, 0.0);
	        NewIterSSolver newSolver = new NewIterSSolver(A, 1.0e-15, 1000000);
	        newSolver.PCGSSOR(RHS, x, 1.5e0, 1);

	        //Arrays.fill(x, 0.0);
	        //newSolver.GMRES(RHS, x, 200);
	        //newSolver.GMRES_SSOR(RHS, x, 1.5, 30, 1);
	        mesh2D.toTecplot("SolutionsT6.dat", x);
//	        System.out.println(MVO.toString(x));
	        
	        double[][] gradX = new double[2][fs.GetNdof()];
	        fs.computeGradient(x, gradX);
	        mesh2D.toTecplot("GradSolutionT6.dat", gradX);
	        System.out.println("Gradient");
	        System.out.println("Dx min:" + MVO.min(gradX[0]) + "\t" + "max: " + MVO.max(gradX[0]));
	        System.out.println("Dy min:" + MVO.min(gradX[1]) + "\t" + "max: " + MVO.max(gradX[1]));
	        
	        double[][] soX = new double[4][fs.GetNdof()];
	        fs.computeSecondDerivative(x, soX);
	        mesh2D.toTecplot("GGradSolutionT6.dat", soX);
	        System.out.println("Gradient2");
	        System.out.println("Dxx min:" + MVO.min(soX[0]) + "\t" + "max: " + MVO.max(soX[0]));
	        System.out.println("Dxy min:" + MVO.min(soX[1]) + "\t" + "max: " + MVO.max(soX[1]));
	        System.out.println("Dyy min:" + MVO.min(soX[2]) + "\t" + "max: " + MVO.max(soX[2]));
//	        ArrayList<Integer>  li = new ArrayList<Integer>();
//	        li.ensureCapacity(100);
//	        for(int i = 0; i < 100; i++) {
//	        	li.add(null);
//	        }
//	        li.add(1, 2);
//	        System.out.println(li.toString());
	        
//	        Mesh2DT3 mesh6 = mg.square2D(2, 2, T3Type.Left);
//	        mesh6.displayBoundarys();
//	        mesh6.extractBoundarysFromElements();
//	        mesh6.displayBoundarys();
		}

	}
